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.SUMMARY 


Although the application of multigrid methods to the equations of 
elasticity and structural mechanics has been suggested, few such applications 
have been reported in the literature. In the present work, multigrid 
techniques are applied to the finite element analysis of the deflections of a 
simply supported Bernoull i-Euler beam, and various aspects of the multigrid 
algorithm are studied and explained in detail. 

In this study, six grid-fineness (or coarseness) levels, with 32 elements 
on the finest grid level, were used to model half the beam. To test the 
multigrid algorithm more severely, random initial approximations were used. 
With linear prolongation and sequential ordering, the multigrid algorithm 
yielded results which were of machine accuracy with work equivalent to 200 
standard Gauss-Seidel iterations on the fine grid. Also with linear 
prolongation and sequential ordering, the V(1,n) cycle with n greater than 2 
yielded better convergence rates than the V(n,1) cycle. 

Derivation of the restriction and prolongation operators was based on 
energy principles. Conserving energy during the inter-grid transfers required 
that the prolongation operator be the transpose of the restriction operator. 
Maintaining an energy balance in the inter-grid transfers also led to improved 
convergence rates. With energy-conserving prolongation and sequential 
ordering, the multigrid algorithm yielded results of machine accuracy with a 
work equivalent to 45 standard Gauss-Seidel iterations on the fine grid. The 
red-black ordering of the relaxation with either linear or energy-conserving 
prolongation yielded solutions of machine accuracy in a single V ( 1 , 1 ) cycle, 
which required work equivalent to about 4 iterations on the finest grid level. 
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INTRODUCTION 


Multigrid methods are known to be efficient solution techniques for 
certain classes of partial differential equations, and their advantages in the 
field of fluid dynamics have been clearly established [1,2]. Although their 
application to the equations of elasticity and structural mechanics has been 
suggested, there have been few attempts to implement multigrid methods in 
these cases. 

Certain basic ideas underlying multigrid methods can be found in the work 
of Southwell [ 3 ]. However, the work of Fedorenko [4] should be regarded as 
the forerunner of the multigrid methods. Fedorenko solved the elliptic 
difference equations, using concepts of error smoothing by Jacobi relaxation 
and calculation of corrections on coarser grids. He also obtained an 
asymptotic estimate of the computational complexity for the Poisson equation 
on a rectangle. Fedorenko’s technique was generalized by Bakhvalov [bj to 
include the advection diffusion equation. At about the same time, Wachpress 
L6 J proposed a two-grid method for elliptic systems. However, multigrid 
methods remained virtually unused for almost a decade until the pioneering 
work of Brandt [7], Nicolaides [8J and Hackbush [9J revived them In the mid- 
1970's. Specifically, Brandt’s early work established the actual efficiency 
of multigrid methods, and demonstrated the capability of these methods to 
treat nonlinear equations, general domains, local grid refinements, and 
solution adaptivity. He demonstrated that the maximum number of computer 
operations required to solve a discrete system of n equations in a Poisson 
problem was on the order of n computer operations. Brandt also showed how a 
local Fourier analysis could be used for the theoretical investigation of 
smoothing rates and optimization of procedures. The work of Nicolaides [8] is 
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the first systematic study of multigrid procedures relating to the finite 
element discretization of elliptic equations. A complete historical 
background of multigrid methods may be found in Stuben and Trottenburg [10j, 
and Hackbush L 11 j . 

Initially, multigrid methods were successfully applied to elliptic 
equations to which they are particularly well suited. Recently, their use has 
been extended to parabolic [12] and hyperbolic [2] equations as well. 
Developments in multigrid methods since the mid-1970's can be found in 
proceedings edited by Hackbush and Trottenburg [13J. Paddon and Holstein [14], 
and McCormick and Trottenburg [15]. 

Many problems in elasticity and structural mechanics are governed by 
elliptic equations; hence, it is quite surprising that their use in these 
fields remains largely unexplored. The present work is an attempt to 
introduce these methods to elastic ians and computational structural 
mechanicists who have had little or no exposure to multigrid techniques. The 
focus of the present work is on a simple, one-dimensional example of a 
Bernoulli-Euler beam for the following reasons: it is amenable to local 

Fourier analysis, yielding smoothing and convergence factors which give an 
idea of the efficiency to strive for; it permits the derivation and 
description of various multigrid elements without the complications of higher 
dimensions; and it serves as a building block for implementing the algorithm 
in higher dimensions. 

Specifically, thi3 paper describes the finite element discretization of 
the relevant beam equations and discusses the performance of the Gauss-Seidel 
iterative technique for the solution of the resulting linear system of 
equations. A particular geometric multigrid scheme for accelerating the 
convergence (usually called the correction scheme) is presented. A general 
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procedure for defining the f ine-to-coarse transfer of residuals (referred to 
as restriction) and coarse- to-f ine transfer of corrections (referred to as 
prolongation) is described from a physical view-point. A discussion of 
storage and work requirements is also included. 

The convergence characterist ics of the solution and the number of 
floating point operations (work) necessary to achieve convergence are 
presented. Results are given for two different orderings of the iterations 
and for two different prolongation schemes. The effect of varying the number 
of iterations at each grid level is also studied. 

NOMENCLATURE 
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Fourier series coefficients, (j =* 1 to N) 
element length on level k, (k * 1 to M) 
semi-bandwidth of stiffness matrix 
flexural rigidity 

value of direct solution at node j, (j = 1 to N g + 1) 
error in iterative approximat ion 
accuracy of the machine 

error in w and 0 solutions at node i, (i * 1 to N + 1 ) 

e 

error normalized by maximum displacement 

force vector 

th ^ „ 

m component of F 

element length on finest level 

structural stiffness matrix 

element stiffness matrix 
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entry in m column, j row of matrix K 
half-length of beam 
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w . 0 . 

J J 


l r l I 


u . 
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total number of grid levels 
concentrated applied moment 
nodal moments 

number of degrees of freedom 

number of elements 

shape functions, (i = 1 to 4) 

number of iterations performed on the downward or first 
leg of the V-cycle 

number of iterations performed on the upward or second 
leg of the V-cycle 

prolongation matrix 

concentrated applied force 

maximum value of q(x) 

loading function applied to beam, (0 £ x £ L) 
restriction matrix 
nodal forces 

i th component of the residual vector 
w and 0 residuals at node j, (j = 1 to N g + 1) 

L 2 ~norm of the residuals 

solution vector 

i til component of U 

approximation to solution vector U 

i tl1 component of V 

work done by concentrated moment 

work done by nodal forces 

work done by concentrated force 

displacement 

arbitrary constants in displacement function 
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0 


rotation 


it strain energy 

Superscr ipts : 

f 0 

, fine and coarse grid functions, respectively 

iteration cycle 

T 

transpose of a matrix or vector 
FORMULATION OF THE PROBLEM 
One-Dimensional Beam Problem 

Consider a simply supported beam of length 2L subjected to the 
distributed loading q(x), as shown in Figure 1. Since the problem is 
symmetric about the center of the beam, only one half of the beam needs to be 
considered. Moment equilibrium requires that 

2 2 

El — j = q Q Lx { L-^L—) 0<x<L (1) 

dx^ 

where El is the flexural rigidity of the beam and q^ is the maximum value of 
the loading function at the center of the beam. 

The boundary conditions for the simply supported beam are w(0) = w(2L) * 0. 

The exact solution for the loading shown in Figure 1 is 

q xL 3 4 2 

W(X) = I20EI - 10 (g-) + 25 } 0 £ x £ L (2) 

To solve the problem using the finite element method, the structure is 
idealized by a number of beam elements. For the two-noded beam element, the 
displacement w and rotation dw/dx (= 6) are used as the nodal parameters. The 


6 



relevant boundary conditions are w(0) = 0 and 6(0 = 0. The element stiffness 


matrix [kj, obtained 

from the 

shape 

functions 

given in 

Appendix A, 

for a beam 

element of length b 

is given 

below. 
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-6/b 2 
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2/b 

-6/b 2 
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The element stiffness matrices are assembled into a structural stiffness 
matrix K and, for any given loading and boundary conditions, the system to be 
solved can be expressed in matrix form as the following equation: 

KU = F (4) 

where K is a positive-definite, symmetric, square matrix, F is the load 
vector, and U is the vector of unknown nodal displacements and rotations. 

Conventional finite element analyses typically use direct solution 
techniques, such as Gaussian elimination or Cholesky decomposition, to solve 
equation (4). With increasing discretization, solution time can become 
critical; thus, an attractive alternative to a direct solution is an efficient 
iterative technique. 

Iterative Solution Techniques 

Two widely used iterative solution techniques are Jacobi iteration and 
Gauss-Seidel iteration. The Gauss-Seidel method is used in the present work. 
Gauss-Seidel Iteration 

Starting with (J 1 as the initial approximation to the vector of unknowns 

s+1 

U, Gauss-Seidel iteration uses the following algorithm to calculate u m , the 
m tl1 component of U: 
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mj J 


(m = 1 ,2 N) (b) 


S th 

where and f are the m components of U and F, and s indicates the 
iteration number. The iteration is continued until the required number of 
iterations has been completed or until the change in the current estimate of 
the solution vector is less than some given tolerance. 

Convergence Study 

To study the convergence of the Gauss-Seidel method, the beam problem 
shown in Figure 1 , with one half the beam modeled with 32 elements, was 
considered. The convergence of the Gauss-Seidel iteration was monitored by 
calculating the L^-norm of the residuals, denoted as |jr||^- In discrete 
form, assuming all elements are of equal length b, the L^-norm can be 
expressed as follows: 



where the residual r. 

i 


is defined as follows: 


r . 
l 


- f. 


N 

l 

j=i 


k . .u . 
ij J 


(i = 1 to N) 


( 6 ) 


( 7 ) 


In the present finite element context, the residuals in each equation can 
be thought of as the unbalanced forces and moments. So that all the residuals 
will have the same dimensions, the moment residuals were normalized by b, the 
element length. A non-dimensional L 2 -norm can then be computed as follows: 


8 



/(q 0 L) (8) 

where r y and r Q are the w and 0 residuals. The L 2 -norm has also been divided 
by the factor q^L (the total load on the beam) to produce a dimensionless 
quantity. 

To start the Gauss-Seidel iteration, a random initial approximation for U 
was chosen. Figure 2 shows the L^-norm plotted versus the number of 
iterations for two procedures. In the first procedure, the iterations were 
performed sequentially from node 1 (at the support) to node 33 (at the center 
of the beam). These results are shown by the solid line in Figure 2. In the 
other procedure, the so-called red-black ordering was used. In this ordering, 
all the even-numbered nodes are designated red and all the odd-numbered nodes 
are designated black. Iterations are performed on all the red nodes first, 
then the black nodes. The results of the red-black ordering are shown by the 
dashed line in Figure 2. For both methods, the L 2 ~norm drops rapidly in the 
first 50 iteration cycles; thereafter, the convergence slows, although the 
red-black ordering produces slightly better results than the sequential 
ordering for large numbers of iterations. However, even after 1600 
iterations, the L^-norm for both is still quite high, indicating that the 
iterative solution is still grossly in error. Appendix B further examines the 
convergence of this problem. 

This behavior is inherent in the Gauss-Seidel method, since Gauss-Seidel 
efficiently removes the high frequency errors but not the low frequency 
errors. However, the low frequency errors, since they are smooth, can be 
approximated on coarser levels where they become higher frequency errors. On 
coarser levels, these high frequency errors can then be removed more 
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efficiently. Even if the errors are still smooth on coarser levels, 
iterations on the coarser levels require much less work due to the decrease in 
the number of unknowns. These ideas are exploited in the multigrid procedure 
in a recursive manner to improve the convergence rate. 

MULTIGRID PROCEDURES 

In this section, the multigrid method is described. In multigrid 
terminology, relaxation is used to mean iterations of some smoothing 
technique, such as the Gauss-Seidel described above; thus, to perform 
relaxations on a set of equations means to iterate using a smoothing method. 
Often, the terms relaxation and iteration are used interchangeably in 
multigrid literature. In the remainder of this paper, the term relaxation is 
used to mean Gauss-Seidel iteration, and n relaxations means n iterations on 
the system of equations. 

The multigrid techniques are iterative methods that use a sequence of 
grids and a simple relaxation scheme such as Gauss-Seidel. As previously 
stated, such relaxation is very efficient for removing the high frequency 
errors on a given grid. The remaining low frequency errors become higher 
frequency errors on coarser grids, where they can be effectively smoothed. 

The key elements of the multigrid procedure are the relaxation technique and 
the coarse grid correction. The multigrid algorithm described here is 
generally called a geometric multigrid method. 

Coarse Grid Correction 

Consider the interplay between a coarse grid and a fine grid, where the 
coarse grid has half as many node points as the fine grid. Let the following 
equation be the discrete finite element representat ion of the given problem on 
a fine grid. 
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( 9 ) 


f f f 

K U = F 


If V f is the approximation obtained on the fine grid after a few relaxation 


f f f 

sweeps, then the error e , = U , - V satisfies the following equation: 

„f f f 
K e = r 

_ f f r r 

where r , the residual on the fine grid, is defined as r = F - K V . 
f f 

If e is a smooth function, e can be approximated by a coarse grid 
function e° which satisfies the following equation: 

K c e c = r G 

r C = R (r f ) = R (F f - K f V f ) 


( 10 ) 


(ID 


where 

Here R is the fine to coarse grid restriction operator. Since the coarse 
grid has half as many grid points as the fine grid, it is more economical to 
solve equation (11) than (10). After solving equation (11) approximately by 
some method, the approximate error e G is used to accelerate convergence of the 


fine grid solution using the following: 
V f «- V f + P(e C ) 


( 12 ) 


Here P is the prolongation operator that interpolates from the coarse grid to 
the fine grid and the arrow «- means that the left side is replaced by the 
right side. 

The process of calculating r C , solving for e°, and interpolating the 
results to the fine grid is called the coarse grid correction. To solve 
equation (9) efficiently, the above procedure is applied recursively using 
multiple grids. This recursive procedure is described in the following 
section. 

Multigrid Cycle 

Figure 3 shows the sequence of grid levels, k - 1 to M, for a simply 

k-1 

supported beam where the element size b, satisfies the relation b^ = 2 h and 
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h is the element size on the finest grid level (k =1). The equation on level 


k is written as follows: 
KV - F k 


(13) 


For k * 1 (the finest level), K 1 * K and F? = F , respectively the given 
stiffness matrix and right-hand side of the original problem. The basic 
algorithm must determine the solution of the algebraic system in equation ( 13 ) 
on the finest grid, k = 1. The fundamental idea is to use the solutions from 
the progressively coarser grid levels to construct other finite element 
solutions in the form of equation 0 3). This procedure is called smoothing: 
in practice, it is carried out by means of relaxation. A commonly used 
technique is the Gauss-Seidel method. 

Starting on the finest level, k - 1, let be an approximation to in 


equation (13), and define the residual on the finest level as r 


F 1 - KV. 


If e 1 denotes the error U 1 - V 1 , the error and the residual are related 

through the following equation: 

„1 1 1 
K e = r 


( 14 ) 


For all subsequent levels (k > 1), the residual r is the restriction of the 


residual from the previous level k - 1 defined by the following equation: 

,k-1 k-1 


k „ , k-1 
r = R (r 


K 


) 


and, therefore, the equation to be solved is now the error equation: 

„k k k 
K e = r 


( 15 ) 


( 16 ) 


Given an approximate solution V to the equation (13), the multigrid 
cycle for improving this approximation can be executed recursively as follows. 

If k = 1 , perform n 1 relaxation sweeps on equation (13), thereby 
obtaining an approximate solution V 1 . Calculate the residual r 1 - F 1 - kV. 
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Down-Sweep : 

For 1 < k < M, repeat steps (1) and (2); for k = M, go to step (3), with 
these steps defined as follows: 

(1) On the current level k, start with an initial approximation 

k 

e =0, and perform n^ relaxation sweeps on equation (To), 
obtaining a new approximation e . 

(2) Restrict the residual to the next coarser level, k + 1, where 

k+1 , k „k k, 

r ■ R(r - K e ), and set k = k + 1. 

(3) If k » M, solve equation (16) exactly, set k ■ M - 1, and start 
the up-sweep. 

Up-Sweep : 

For M > k > 1, repeat steps (4) and (5); when k = 1, go to step (6), 
where these steps are defined as follows: 

k 

(4) Update the approximation to the error e on the current level 

by adding the prolongated correction of the error from the 

previous coarser level: 
k k k+1. 

e <- e + P ( e ) ( 1 '/ ) 

(5) Perform n^ relaxations on equation (16), starting with the 
updated error e from equation (17). This will result in an 

k 

improved approximation to e . Set k = k - 1 . 

(6) For k = 1 , calculate an updated approximation to the solution 

1 1 2 

V <- V + P(e ), and perform n^ relaxations on equation (14), 
obtaining an improved approximation to V 1 . 

The cycle described above in steps (1) through (6) is called a V-cycle, 
indicated by the notation Vtn^.n^), where n^ indicates the number of 
relaxation sweeps performed at each level on the first or downward leg of the 
V-cycle and n^ indicates the number of relaxation sweeps performed at each 
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level on the second or upward leg of the cycle. Figure 4 shows a schematic 


drawing of a V-cycle for six grid levels. A flowchart for one V-cycle is 
shown in Figure 5. 

Restriction and Prolongation Matrices 
From the foregoing description of the correction scheme, it is apparent 
that the keys to the multigrid procedure are the restriction and prolongation 
operators. These operators are related to each other as explained below. 

From equations (16), the residual relations on the coarse and fine grid 


are written as follows: 


„c c c 
K e = r 


„f f f 
K e = r 


( 18 ) 


If e is smooth, it can be represented on a coarser level. In that case, 
the coarse and fine grid errors can be related using the prolongation operator 
P in the following equation: 


e f - Pe° 


( 19 ) 


f c 

and the residuals r and r can be related through the restriction operator R 


as follows: 


r G = Rr f 


( 20 ) 


In equations (18), the errors, e and e , can be thought of as representing 
the displacements of the respective grids, while the residuals, r G and r f , can 
be thought of as the forces required to produce these displacements. Thus, 
the strain energy of the coarse and fine grids can be expressed in terms of 
the errof equations by the following expressions: 


1 , c , T c 

2 (e ) r 


1 , f.T f 

2 (e ) r 


(21 ) 


From equations (18), (19) and (20), the strain energies of the coarse and 
fine grids may be rewritten as follows: 


1*4 



( 22 ) 


, c . T_ 
(e ) R 


f 

TT 


(e c )V 


f 

r 


c 

If the strain energy is to remain constant during intergrid transfers, it 

f c f T 

must equal ir . Thus, for arbitrary values of e and r , R must equal P to 

maintain an energy equivalence between grid levels. 

One-Dimensional Beam Example 

Before applying the multigrid cycle to the solution of the beam problem, 
the restriction matrix R and the prolongation matrix P must be defined. 
Restriction Matrix 

The restriction matrix R is defined to transfer the loads and moments 
from a fine grid to a coarse grid. To form this matrix, the consistent load 
relations must be formulated for transferring the loads and the moments. With 
equations (A2) and (A3), it is possible to determine the nodal forces in the 
beam element equivalent to a concentrated force and moment applied at the 
center of the element. The details of this formulation are given below. 

Concentrated force . The nodal forces equivalent to a concentrated force 
Pq applied at the center of an element of length b are found by equating the 
work done by the concentrated force to the work done by the unknown nodal 
forces. 

Wp, the work done by the concentrated force, is equal to the force 
times the displacement of that force: 

Wp = Pq x w(x=b/2) (23) 

The substitution of x * b/2 into equations (A2) and (A3) yields the following 
result: 

w p = p 0 x( W i /2 + e^b/8) + Wp /2 - e 2 (b/8)) (24) 
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Similarly, the work done by the nodal forces acting at the two end nodes 


can be written as follows: 


“n ‘ Vl * * M 1 8 1 4 M 2 6 2 

where are the nodal forces and are the nodal moments. 

Wp should equal for any arbitrary values of w 1 , w 2 , 0 1 , and 0^. 
the consistent forces are 


(25) 
Thus , 


R 1 - R 2 - V 2 


M 1 =P Q b/8 


M = -M 
2 1 


(26) 


Concentrated moment . As before, the nodal forces equivalent to a 
concentrated moment M Q applied at the center of the beam element are found by 
equating the work done by the concentrated moment to the work done by the 
unknown nodal forces. 

W M , the work done by the concentrated moment, is equal to the product of 
the moment and the rotation caused by that moment: 

W M = M 0 x 0(x=b/2) (27) 


To determine 0 in an element, equation (A2) must be differentiated with 


respect to x, yielding the following equation: 

e = dw = L _ 6x + 6xf + C1 _ 3x1-. 

dx 1 L 2 k 3 1 b 2 J 

b b b 


(28) 


+ w L 6 x _ 6x 2 j + L _ 2x + 3x 2 

2 L b 2 b 3 J V b b 2 J 


0 < x S b 


Substituting x = b/2 into equation (28) and simplifying, the work due to the 
concentrated moment M can be written as follows: 

W M = M 0 x[( " 3/2b)w i ' V 4 + (3/2b)w 2 ' 6 2 / ' 4j (29) 
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As for the concentrated force, the work due to the nodal forces Is 
written as follows: 

W n = R 1 w 1 + R 2 w 2 + M 1 6 1 + M 2 6 2 (30) 

where , R 2 are the nodal forces and M 1 , M 2 are the nodal moments. 

Again, W M should equal for any arbitrary values of w 1 , w 2> 9^ and 0 . 
Thus, the equivalent nodal forces for the case of the concentrated moment can 
be expressed as follows: 


R 1 - - 3M Q /(2b) 


R 2 ' ~ R 1 


M = - M /4 


M = M 
2 1 


(3D 


Consistent residual transfer . The consistent nodal forces and moments 
given in equations (26) and (31) are used to determine the residual weights 
fof use in the energy-conserving inter-grid transfer. 

To form the restriction matrix using the force and moment weights, 
consider two grid levels as shown in Figure 6. Here the fine grid has three 
node points and two elements, while the coarse grid has two node points and 
one element. The lengths of the elements in the fine and coarse grids are b/2 
and b, respectively. The restriction matrix R is defined to restrict the 
forces and moments from the fine grid to the coarse grid. The transfer of the 
forces and moments is shown in Figure 6 and defined below. 


r G = R[r f ] 


where 


r r - IK; Mj Rj Mj R k M k ) T r° - IR,, R k , M,,) 7 
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r 1 

0 

1 /2 

-3/ (Mb) 

R = 

0 

1 

b/4 

-1/4 


0 

0 

1/2 

3/ (4b) 


- 0 

0 

-b/4 

-1/4 


0 0 “ 

0 0 

1 0 

0 1 - 


( 32 ) 


This restriction matrix is valid only when the element size is doubled at each 
coarser grid level. 

Prolongation Matrix 

The prolongation matrix defines the displacements on the fine grid given 
the displacements on the coarse grid. To obtain this matrix, again consider 
the two grid levels shown in Figure 6. Here the displacements w and 8 are 
known on the coarse grid at nodes i T and k 1 . From these displacements, the 
displacements at nodes i, j and k on the fine grid can be found in two ways. 

In either method, the displacements at nodes i and k of the fine grid are the 
same as the displacements of nodes i f and k 1 on the coarse grid. The 
displacements at node j on the fine mesh are found in two ways. One method 
assumes that the displacements at node j are the average of the displacements 
at nodes i T and k T on the coarse mesh. This amounts to linear interpolat ion 
as defined below: 


w C « PLw f J 


where 


w = { w . 


0 . w . 0 . w. 0, ) 

i J J k k 


1 

0 

1/2 

0 

0 

0 


0 

1 

0 

1/2 

0 

0 


0 

0 

1/2 

0 

1 

0 


0 

0 

0 

1/2 

0 

1 



V 


w, 


v 1 


T 


( 33 ) 
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In the other method, the element shape functions given in equations (A2) 
and (A3) are used for the interpolation. The displacements at the center of 
the element in the coarse mesh are found by substituting x = b/2 into 
equations (A2) and (A3). These are the displacements at node j on the fine 
mesh. The prolongation matrix thus obtained is given below: 



- 1 

0 

0 

o - 


0 

1 

0 

0 

p - 

1 /2 

b/4 

1 /2 

-b/4 


-3/(4b) 

1 /4 

3/ (4b) 

-1/4 


0 

0 

1 

0 


0 

0 

0 

1 J 


The prolongation matrix obtained using the element shape functions is 
exactly the transpose of the restriction matrix given in equation (32). Thus, 
the use of this prolongation matrix maintains an energy equivalence between 
grid levels. Again, this prolongation matrix is valid only when the element 
size is doubled at each coarser grid level. 

Although the restriction and prolongation operators derived here are for 
a one-dimensional problem, the procedure used here can be applied to two- or 
three-dimensional problems. 

Storage and Work Requirements 

The storage requirements necessary to solve the equation KU=F using 
direct solution techniques are easily calculated from the total number of 
degrees of freedom in the problem and the semi-bandwidth of the stiffness 
matrix K. Since the multigrid method is being proposed as an alternate 
solution technique, the storage requirements of the multigrid method should 
also be addressed. 
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Obviously, the multigrid method requires more storage space than direct 
solution techniques if all entries within the bandwidth are stored. The 
original stiffness matrix must be stored as the finest grid level and 
additional grid levels must be accommodated. If the finest grid has N g 
elements, the coarser levels have N e /2, N g /4, etc., elements. Thus, the 

total number of elements for M levels is 

N + N/2 + N/4 + N/ 8 + ... + N /(2 M_1 ) = 2N (1 - ^ 77 ) < 2N (35) 
e e e e e e 2 M e 

The total number of degrees of freedom is computed in a similar fashion. In 
the present one-dimensional case, for N g elements on the finest grid, there 
are (N g +1) nodes, each with 2 degrees of freedom. Hence, the total storage 
required, assuming a constant semi-bandwidth B, can be found from the 
following expression: 

N N 

2B(N e + 1) ♦ 2B( + D ♦ 2B(jp + 1 ) + ... = 2B(M + 

2 N e {1 ~ ^M }) < 2B(2N e +M) ( 36 ) 

Thus, the maximum storage required is always less than 2B(2N g +M). 

This computation assumes the storage of all entries within the bandwidth, 
which may result in the storage of several zero entries. These zero entries 
must be stored in direct solution techniques since during the solution phase 
some zero entries may become nonzeroes. Examination of the Gauss-Seidel 
iteration in equation (5) shows that only the nonzero terms of the matrices 
need to be stored. However, storing only the nonzero entries requires the 
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additional storage of pointer arrays for each row of the stiffness matrix. 

Many finite-element stiffness matrices are sparsely populated and, in such 
cases, the storage of only the nonzero terms, with the pointer arrays, can 
considerably reduce the storage requirements compared with banded or profile 
storage schemes. 

With the definition of the coarser grids used here, the total storage 
required by the stiffness matrices of all levels should be less than twice 
that needed by the finest level. Additionally, the solution vectors for each 
grid must be stored, as well as the load vectors, or right-hand sides. 
Depending upon the structure of the problem, it may be advantageous to store 
the restriction and prolongation matrices as well. 

The work (number of floating point operations) done in a single Vln^n^) 
cycle is equivalent to 2(n^ + n^) times the work of a single relaxation on the 
finest grid level. This is determined as follows. In the VCn^n^,) cycle, 
n.j+ n^ relaxations are done on the finest level while the total sum of the 
relaxations done on the lower levels is always less than or equal to n^+ n^ . 
Because each subsequently coarser grid level has half as many degrees of 
freedom as the one before, the total number of degrees of freedom in all 
levels is always less than twice the number in the finest grid (see equation 
( 36 )). The work done iterating on each level is directly related to the 
number of degrees of freedom that must be relaxed. Thus, the total work done 
in a single V(n 1 ,n 2 ) cycle is less than the work of 2(n 1 +n^) relaxations on 
the finest grid level. 

RESULTS AND DISCUSSION 

In this section the performance of the multigrid algorithm in solving the 
simply supported beam problem of Figure 1 is evaluated. Six grid levels are 
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used, where the finest mesh (level 1) has 32 elements modeling half the beam, 
the next coarser grid (level 2) has 16 elements modeling half the beam, and so 
on until the coarsest grid (level 6), which has one element modeling half the 
beam. 

The two prolongation operators mentioned earlier are considered: a 

simple linear interpolation and a prolongation operator based on a constant 
energy in each grid level. 

Results are presented for a VCn^n^) cycle. The convergence of the 
solution and the work necessary to achieve convergence are presented. 

Convergence of the solution was monitored with the L 2 -norm of the 
residuals, defined in equation (8). While the displacements converge to the 
exact solution within machine accuracy, the L^-norm, because of its 
definition, converges with less accuracy than the displacements. Appendix C 
further explains this behavior. 

Results are given first for sequential ordering with linear prolongation, 
for varying values of n^ and and then with the energy-conserving 
prolongat ion. Finally, the effects of red-black ordering are presented. 

Sequential Ordering 

Linear Prolongation 

Figure 7 presents the L 2 -norm plotted against an increasing number of V- 
cycles for both random (solid line) and zero (dashed line) initial approxima- 
tions. These results are for a V(2,1) cycle where the coarsest level was 
solved exactly. As seen in the figure, at the end of 20 V-eycles, the L^-norm 
is on the order of 10 for the random initial approximation and on the order 
of 10 ^ for the zero initial approximation, indicating that the multigrid 
solution is close to the direct solution obtainable on the finest grid level. 
The r.-;ues of convergence (indicated by the slope of the curves) are almost 
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identical for both the zero and random initial approximations. To test the 
algorithm more severely, only the random initial approximation is used in the 

remainder of this work. The work necessary to achieve the accuracy shown in 

Figure 7 is 2(n^+ n^) times 20 V-cycles, or equivalent to the work of 120 
relaxations on the finest grid level. 

Varying n 1 and n 2 

To study the effect of varying the number of relaxations on either leg of 
the V-cycle, two cases with linear prolongation were considered. In the first 
case, the number of relaxations on the downward leg of the cycle was held 
constant at = 1 , while the number of relaxations on the upward leg of the 
cycle was increased from 1 to 5 (n 2 ■ 1 to 5). Figure 8 shows the L^-norm as 

a function of increasing n^. As expected, the rate of convergence improves 

-7 

with increasing n^. For the V(1,5) cycle, the L^-norm is on the order of 10 
after about 18 V-cycles. The work of 18 V(1,5) cycles is equivalent to that 
of 216 relaxations on the finest grid level. 

The second case, n^ = 1 while n 1 was increased from 1 to 5, is shown in 
Figure 9. As in the previous case, the rate of convergence improved with 
increasing numbers of relaxations. For the V(5,1) cycle, the L -norm is on 
the order of 10 ^ after about 18 V-cycles, compared to 10 ^ for the V(1,5) 
cycle of the previous case. As before, the work of 18 V(1,5) cycles is 
equivalent to that of 216 relaxations on the finest grid level. A comparison 
of Figures 8 and 9 shows that the V(1,n) cycle yields better convergence than 
the V(n,1) cycle when n is greater than 2. This can be explained in the 
following manner. 

A prolongation from level k to k-1 introduces high frequency errors at 
level k-1. This is illustrated schematically in Figure 10 for linear 
interpolation. Level k represents a coarse mesh with 2 elements, AC and CL. 
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Level k-1 represents a fine mesh with 4 elements, AB, BC, CD, and DE. Using a 
linear prolongation, the displacements on level k-1 at B and D are found as 
= (w ft + w c )/2 and w D = (w + w £ )/2 with similar expressions for the 
rotations at B and D. The errors in the displacements at B and D are larger 
than at C and E, as shown in the Figure 10(b). The total error can be 
represented as a low frequency (smooth) error plus high frequency 
(oscillatory) errors. These high frequency errors can be easily eliminated on 
level k-1 with a few relaxations. The V(1,n) cycle does precisely this, and 
hence the better convergence compared with the V(n,1) cycle. 

Energy-Conserving Prolongation 

Figure 11 shows the variation of the L^-norm as a function of the number 

of V-cycles for the V(1 ,n 1 ) cycle with n 1 =1,2 5 and P * R T (energy- 

conserving prolongation). The larger the value of n 1 , the faster the solution 
converges. For n 1 equal to 3, 4 or 5, the multigrid results converge to the 
exact solution within the accuracy of the machine in approximately 8 V-cycles 
with a work equivalent to about 96, 80 or 64 relaxations on the finest grid. 
For n 1 equal to 1 and 2, slightly more V-cycles were required for convergence 
for a work expenditure equivalent to 48 and 60 iterations, respectively, on 
the finest grid. Thus, although for n^ equal to 1 or 2, more V-cycles are 
required for convergence, less work is needed than when n 1 is equal to 3, 4 or 
5 . 

T 

The energy-conserving prolongation (P = R ) was also used with the V(n,1) 
cycle. Two cases are compared in Figure 12. The results for the V(n,1) cycle 
are virtually identical to those for the V(1,n) cycle. This is expected since 
the energy-conserving prolongation introduces little or no high frequency 
errors from the coarse level to the fine level. 


24 



Figure 13 compares the exact solution to the multigrid solution for 
energy-conserving prolongation with the V(2,1) cycle. Figure 13 shows e w and 
e , the errors in w and 9, respectively, where e and e n are calculated as 

U Wo 


e 

w 


w(x) 


- w(x) 


w 

max 


exact 


9 ( X ) - 0 ( x ) 
9 

max 


exact 


Here w and 8 are the maximum values of the exact solution. Figure 13 
max max 

shows e and e at the end of the first, third, fifth, and seventh V-cycles. 

W v 

After each V-cycle, the maximum error drops significantly, and after the 
seventh cycle, the multigrid solution has converged to the exact solution 
(within machine accuracy). 

Red-Black Ordering 

To study the performance of the red-black ordering scheme in the V(n 1 , n^) 
multigrid cycle, all the even-numbered nodes were designated as red and all 
the odd-numbered nodes were designated black. At each grid level, all the red 
nodes were relaxed first, followed by the black nodes, in a V(1,1) cycle. For 
both the linear and energy-conserving prolongations, the solution converged 
within machine accuracy in a single V-cycle. Thus, with the red-black 
ordering, convergence is obtained with a work equivalent to only 4 relaxations 
on the fine grid. The excellent performance of the red-black ordering scheme 
for this one-dimensional problem is expected and can be explained by a cyclic 
reduction method. 

In the cyclic reduction method, the displacement and rotation at node i 
(red node) can be expressed in terms of the displacements and rotations at 
nodes i-1 and i+1 (black nodes), and the node i displacement and rotation can 
be eliminated from the system. This can be done for all the red nodes, thus 
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producing a much smaller set of equations- The same procedure can be employed 
repeatedly until only one or two nodes remain. The solution for the original 
system Is then found by reversing the process. This type of process, which 
can be applied only to one-dimensional problems, always produces an exact 
solution in one cycle. For this problem, the current process with red-black 
ordering is exactly equivalent to this cyclic reduction scheme, and thus 
converges in one cycle. 


CONCLUDING REMARKS 

Multigrid procedures were developed for use in one-dimensional finite 
element analysis. Several aspects of the multigrid procedure were studied and 
explained in detail. 

The key elements of the multigrid method are the restriction and the 
prolongation operators. These operators were derived based on energy 
principles. A general procedure for obtaining the restriction matrix was 
outlined. This procedure considered the residuals from the iterative solution 
as the loads on the fine grid and transferred these loads, In an energy- 
conserving manner, onto the coarse grid. The procedure outlined here can be 
used with any multigrid finite element procedure. 

Various aspects of the V-cycle multigrid algorithm were studied In the 
analysis of the deflections of a one-dimensional, simply supported Bernoulli- 
Euler beam. From this study, for six grid levels with 32 elements in the 
finest grid, the following conclusions were drawn. 

1. With linear prolongation and sequential ordering, the multigrid 
algorithm yielded results which were of machine accuracy with work 
equivalent to about 200 standard Gauss-Seidel relaxations on the 
finest grid. 
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2. With linear prolongation and sequential ordering, the V(1,n) cycle 
with n greater than 2 yielded better convergence rates than the 
V(n,1) cycle. 

3. Maintaining an energy balance during the inter-grid transfers 
required that the prolongation operator be the transpose of the 
restriction operator. 

4. The multigrid algorithm showed improved convergence rates when energy 
was conserved in the inter-grid transfers. With the energy- 
conserving prolongation and sequential ordering, the multigrid 
algorithm yielded results which were of machine accuracy with a work 
equivalent to about 50 standard Gauss-Seidel relaxations on the 
finest grid. 

5. The red-black ordering of the relaxation with either the linear or 
energy-conserving prolongation yielded solutions with machine 
accuracy in a single V ( 1 ,1 ) cycle. This is equivalent to performing 
about H relaxations on the finest grid level. 
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APPENDIX A - BEAM ELEMENT SHAPE FUNCTIONS 


For an element with two nodes (four unknowns), the deflected shape w can 

be assumed to be given by a cubic function, i.e., 

2 3 

w = a 1 + a 2 x + a^x + a^x 0 £ x S b (A1) 

Equation (AT) can be written in terms of shape functions (i = 1 to 4) 
corresponding to Wj and 9^ (j =* 1,2) as follows: 

w(x) = w ] N 1 + 6 1 N 2 + w 2 N 3 + 0 2 N^ (A2) 

where 


N 


1 



N 


3 



x - 


2x 


N„ = - z — + 


0 < x £ b 


(A3) 


Here b is the element length. 

The element stiffness matrix, obtained from these shape functions, is 
given in the text (equation (4)). 


30 



APPENDIX B - CONVERGENCE OF THE GAUSS-SEIDEL METHOD 


The Gauss-Seidel iterations efficiently remove the high frequency errors 
but not the low frequency errors. This appendix demonstrates the efficiency 
of Gauss-Seidel iterations in removing the high frequency errors and discusses 
the convergence of the Gauss-Seidel method. 

Figures B1 and B2 show the efficiency of the Gauss-Seidel method in 
removing the high frequency errors. A random initial approximation was used 
for these figures. The displacement w at each node, normalized by the exact 
displacement at the center of the beam, is shown at the end of the 1 st , 10 th , 
t>0 th , 1000 th , 10000 th , and 20000 th iteration in Figure B1 . Similar 
information for the rotation 6 at each node is given in Figure B2. For both 
the displacement w and the rotation 0, the majority of the high frequency 
errors have been smoothed out by the 10 th iteration. The remaining smooth 
errors are not handled efficiently by the Gauss-Seidel method. 

Figures B3 and B4 present the iterative solutions for the displacements 
and rotations for a zero initial approximation. The displacement w at each 
node, normalized by the exact displacement at the center of the beam, is shown 
at the end of the 500 th , 1000 th , 5000 th , 10000 th , 20000 th , and 30000 th 
iteration in Figure B3. Similar information for the rotation 0 at each node 
is given in Figure Bk. Here again the Gauss-Seidel method performs very 
inefficiently. 

In fact, the performance of the Gauss-Seidel method is so poor for this 
example that even after 30000 iterations, the iterative solutions are still 
grossly in error. Figures B5 and B6 compare the exact solution to the 
iterative solution for two different initial approximations. Figure BS shows 
w, the displacements, and Figure B6 shows 6, the rotations after 30000 
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iterations. In both figures, the solid line represents the exact solution 
while the dashed and dash-dot line represent the zero and random initial 
approximations, respectively. Even after 30000 iterations, the magnitudes of 
both iterative solutions are grossly in error and, obviously, the nature of 
the initial approximation still has a large effect on the iterative solution. 

The very slow convergence in this problem is also due, in part, to the 
nature of the stiffness matrix for the beam example. The Gauss-Seidel 
iterative method works best for diagonally dominant matrices. However, as 
seen in equation (3)» the beam stiffness matrix contains off-diagonals terms 
of the same or higher order than the diagonals terms. Hence, the poor 
performance of the Gauss-Seidel in this example is not unexpected. 
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APPENDIX C - ERROR ANALYSIS 


This appendix examines the relationship between the error in the 

displacements and the error in the L^-norm of the residuals. As stated 

previously, although the displacements appeared to have converged to the 

-1 4 

direct solution within machine accuracy (10 ), the L^-norm of the residuals 

-9 

was on the order of 10 . Figure Cl shows the errors in w and 0b, normalized 

by the maximum displacement, where the error in the displacement is computed 

as the difference between the iterative solution displacement and the direct 

solution displacement. Figure C2 presents the normalized force and moment 

residuals along the length of the beam. As seen in Figure Cl , the 

-9 

displacement errors are smooth and on the order of 1 0 , and while the 

-9 

residuals (Figure C2) are also on the order of 10 , they are oscillatory. 

Consider the equation for the i tia residual: 


r . 

l 



N 

l 

j-i 


k . .v. 
ij J 


(Cl ) 


where is the current approximate solution. The term can be expressed as 


follows: 


v . - d - e . 
J j J 


(C2) 


where d is the direct displacement solution at node j and e. represents the 
J . J 

error at node j . 


Since by definition f « E k d., equation (Cl) can be rewritten as 

i j-1 ij J 


follows: 
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I 
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k . .e . 
ij J 


(C3) 
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or, 


r, = )" k . ,e , £ > Ik, .lie. I 
1 j=i J J ^ 1 1 iJ 1 1 J 1 


(C4) 


If all the errors e^ (j =1, 2, N) are on the order of e m> the 

accuracy of the machine, then the following is true: 


' S e y k . . 
i m ; L , ij 


(C5 ) 


j = 1 


To keep the dimensions the same, in the current problem of a one- 
dimensional beam, the error in w is denoted by e while the error in 6 is 

w 

denoted by be . Equation (C5) then reduces to the following: 


r < e ^ [12 + 6 + 24 + 12 + 6] 
»i m b 3 


(C6) 


r /b S e — L6 + 2 + 8 + 6 + 2] 

8 i m b 3 

So that the residuals of the forces and moments will have the same dimensions, 
the moment residuals have been divided by b, the element length. 

Equations (C6) can be written in a nondimensional form as follows: 


r / ( q L) £ e 8(L/b) 3 
w . u 
1 

(C7) 

r /(q bL) S e 24(L/b) 3 
e . u 
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where e = and d is the maximum displacement from the direct 

m max max 

solution of the stiffness equations. 

-1 4 

Thus if e m is on the order of the machine accuracy (e m - 10 ), then the 

nondimens ional residuals for the w and 9 degrees of freedom for a beam with 32 

_8 

elements can be at most 1.96x10 and 0.786x10 , respectively. When the 

residuals are smaller than these maximums, as in Figure C2, the solution 
displacements have converged within machine accuracy, and no further 
improvements can be expected. 
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simply supported beam 



beam element 


i q 0 x/L 0 ^ x < L 

q(x) = 

I q 0 (2*x/L) L < x < 2L 
distributed loading 

Figure 1 . Simply supported beam under distributed loading. 
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for sequential and red-black ordering of. Gauss-Seldel 




simply supported beam 
Grid Element 



3 L/8 
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5 L/2 • 

6 L «— — — — • 

Figure 3. Sequence of grids for simply supported beam. 
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Figure 4. V-cycles for six levels, V(n 1t n 2 ). 
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Figure 5. Flowchart of amltlgrid cycle. 








Rj. = Rj + Rj/2 - 3Mj/(4b) R k -= R k + Rj/2 + Mj/(4b) 
Mj- = Mi + Rjb/4 - Mj/4 M k - = M fc - Rjb/4 - Mj/4 


Figure 6. Consistent transfer of force and moment residuals. 
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Error in prolongated solution at level k—1 (fine). 


Figure 10. High frequency errors due to linear prolongation. 
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(a) Error after one V-cycle. 


Figure 13 


e. 



Error in multigrid solution for energy-conserving prolongation and 
V (2 , 1 ) cycles. 
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Figure 13. Continued. 
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Figure B1 . Exact and iterative solutions for deflections with random initial 
approximation. 
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Figure B2. Exact and iterative solutions for rotations with random initial 
approximation. 
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Figure BM. Iterative solutions for rotations with zero initial approximation. 
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Figure B5. Exact and iterative solutions for deflections after 30000 
iterations. 
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Figure B6. Exact and iterative solutions for rotations after 30000 
iterations. 
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